Asian Herpetological Research 2012, 3(2): 103-113 
DOI: 10.3724/SP.J.1245.2012.00103 


Testing Hypotheses of Pleistocene Population History Using 
Coalescent Simulations: Refugial Isolation and Secondary 
Contact in Pseudepidalea raddei (Amphibia: Bufonidae) 


Bingjun DONG"’’, Jing CHE’, Li DING’, Song HUANG’, Robert W. MURPHY” ®, Ermi 
ZHAO'* and Yaping ZHANG” 


' Key Laboratory for Bio-resources and Eco-environment (Ministry of Education), College of Life Sciences, Sichuan 
University, Chengdu 610064, Sichuan, China 

? College of Chemistry and Life Sciences, Shenyang Normal University, Shenyang 110034, Liaoning, China 

* State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of 
Sciences, Kunming 650223, Yunnan, China 

* Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu 610064, Sichuan, China 

` Department of Biology, Huangshan University, Huangshan 245021, Anhui, China 

° Centre for Biodiversity and Conservation Biology, Department of Natural History, Royal Ontario Museum, 100 
Queen s Park, Toronto, ON, M5S 2C6, Canada 


Abstract The impact of the Quaternary glaciation on eastern China’s local fauna and flora is a topic of considerable 
interest. We use mitochondrial DNA (mtDNA) sequences and coalescent simulations to test two general biogeographic 
hypotheses related to the effect of the Pleistocene climatic fluctuations on a widespread, eastern Chinese amphibian, 
Pseudepidalea raddei. Genealogical reconstructions are made and they detect major western and eastern lineages, 
which overlap in northwestern China, and possibly indicate the secondary contact of the populations that had entered 
the region from separate glacial refugia. Coalescent tests rejected alternative hypotheses of fragmentation of either a 
widespread ancestor or panmixia. The tests instead supported the hypothesis of geographic isolation and a remarkable 
dispersal pattern in one of the lineages. Though the Pleistocene climatic events are known to have affected the historical 
distributions and intra-specific divergence of Chinese squamates, coalescent and non-coalescent demographic analyses 
indicated that the toad P. raddei was not adversely affected by glacial cycling. Presumably, an increase in the amount 
of climatically mild habitats in East Asia is due to the development of monsoons since the Mid-late Pleistocene is 
responsible for the relatively mild effects. 


Keywords Mongolian toad, Pleistocene refugia, phylogeography, biogeography 


1. Introduction the microevolution of East Asian species has focused 
largely on the locations of refugia during the Pleistocene 
Quaternary climatic oscillations, punctuated by the 
Pleistocene glaciations, caused massive changes to 
the distribution of species in the Palaearctic realm 


(Hewitt, 2000; Schmitt, 2007). Recent research on 


glaciations and dispersal routes during the postglacial 
period (Stéck et al., 2006; Zhang et al., 2008; Song 
et al., 2009; Ding et al., 2011). Much of the research 
has centered on tropical and subtropical species with 


geographic distributions that sometimes include the 
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Palaearctic (Ding et al., 2011). Regardless, few detailed 
analyses of the populations from the Palearctic are 
available (Bauert et al., 1998). 

Widely distributed amphibians, being terrestrial 
poikilotherms with narrow habitat requirements and 
limited dispersal potentials, are useful model animals for 
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population genetic studies. Amphibians are frequently 
used in ecological and phylogenetic studies, yet no 
intraspecific studies of widespread East Asian Palearctic 
amphibians exist (Hu et al., 2007). 

The Mongolian toad was classified to Bufo (Zhao 
and Adler, 1993), but Frost et al. (2006) considered that 
it should be provisionally assigned to Pseudepidalea. 
Studies on the Mongolian toad focus on ecological (Zhou 
and Song, 1997; Yu and Zhao, 2007), karyotypic and 
genetic analyses (Shang et al., 1982; Yang, 1982; Wang 
et al., 1997; Dong et al., 2004). However, none of these 
studies examine population structure using molecular 
methods. 

Herein, we use the Mongolian toad, Pseudepidalea 
raddei, to investigate the historical demographic 
and genetic effects of recent ice ages on East Asian 
amphibians. This toad occurs from eastern China and the 
Korean Peninsula westward to Mongolia and southern 
Siberia in Russia. The northernmost distribution reaches 
Lake Baikal and the Amur River (Heilongjiang) basin. 
The greatest extent of the species’ distribution occurs in 
China where it ranges throughout much of the north. The 
southernmost extent for a population occurs in Jiangsu 
and Anhui, China (Zhao and Adler, 1993). Because 
populations of this species occur between northern and 
southern borders of the north temperate zone, the species 
is an ideal model animal for use in testing alternative 
hypotheses regarding the persistence of animals in the 
palearctic realm. We gathered and analyzed partial 
sequences of the highly variable mtDNA control region, 
from across the range of P. raddei in China (Figure 1), 
to obtain insights into the population genetic structure of 


the species. We also used coalescent simulations to test 
whether the populations in different areas are descended 
from a wide-ranging common ancestor whose range 
became fragmented (rang fragmentation hypotheses) or 
alternatively, if the pattern of diversification is consistent 
with the biogeographic hypotheses (geology/refugia 
hypotheses) that intraspecific divergences are related 
to geographic distribution and that climatic oscillations 
during the Pleistocene promoted the split of populations 
of this species. 


2. Materials and Methods 


2.1 Sampling A total of 185 specimens of P. raddei 
were collected from 19 localities including most of the 
distributional range of this species in China (Figure 1; 
Table 1), but it was not sampled in Mongolia, Russia, 
and Korea. Voucher specimens were deposited in the 
Kunming Institute of Zoology. Two other toads, P. 
pewzowi and Bufo japonicus (retrieved from GenBank) 
were used as outgroup taxa. Specimens were collected by 
hand and killed following animal use protocols approved 
by the Kunming Institute of Zoology Animal Care 
Committee. 


2.2 Laboratory protocols Total DNA from either 
muscle or liver tissue was extracted according to the 
standard phenol—chloroform method (Sambrook et al., 
1989). Then, DNA was dissolved in ddH,O and stored 
at -20 °C. Primers of the mtDNA control region for P. 
raddei were designed according to the partial sequences 
of the control region from Goebel et al. (1999). Primers 


Table 1 Collection data for specimens of Pseudepidalea raddei included in this study in China. 


Site Locality Voucher No. No. Coordinates Altitude (m) GenBank Accession No. 
1 Harbin, Heilongjiang (HB) YPX838-856 19 N45.41°, E126.15° 201 JQ364649-JQ364667 
2 Xiaoling, Heilongjiang (XL) YPX858-859 2 N45.20°, E127.18° 254 JQ364717—-JQ364718 
3 Shenyang, Liaoning (SY) YPX861-874 14 N41.50°, E123.25° 47 JQ364668-JQ36468 1 
4 Balinzuoqi, Neimeng (BL) YPX88 1-895 15 N43.58°, E119.22° 524 JQ364634-JQ364648 
5 Jining, Neimeng (JN) YPX901-917 17 N41.01°, E113.05° 1447 JQ364700-JQ364716 
6 Hulunbeier, Neimeng (HL) YPX920-937 18 N48.19°, E117.34° 520 JQ364682—JQ364699 
7 Wulashan, Neimeng (WS) YPX02284—02290 7 N41.51°, E109.02° 1428 JQ364587—-JQ364593 
8 Daqindao, Shandong (DD) YPX970-986 17 N37.58°, E120.38° 56 JQ364552-JQ364568 
9 Liulin, Shanxi (LL) YPX939-946 8 N37.25°, E110.53° 783 JQ364624—-JQ36463 1 
0 Binxian, Shaanxi (BX) YPX95 1-966 16 N35.02°, E108.04° 843 JQ364569-JQ364584 
1 Yulin, Shaanxi (YL) YPX949-950 2 N38.17°, E109.44° 1065 JQ364632-JQ364633 
2 Huimeng, Henan (HM) YPX991—1000 10 N34.48°, E112.37° 123 JQ364542-JQ364551 
13 Jiayuguan, Gansu (JG) GS200-209 10 N39.38°, E98.21° 1663 JQ364597—-JQ364606 
14 Jingyuan, Gansu (JY) YPX1004—1016 13 N36.33°, E104.41° 1450 JQ364607—-JQ364619 
5 Yuzhong, Gansu (YZ) GS192-195 4 N35.50°, E104.06° 1978 JQ364620—-JQ364623 
6 Guyuan, Ningxia (GY) NX41-43 3 N36.02°, E106.14° 1782 JQ364594-JQ364596 
7 Helanshan, Ningxia (HS) YPX1002—1003 2 N38.43°, E106.27° 1150 JQ364585—JQ364586 
8 Hainanzhou, Qinghai (HZ) YPX02152-02154 3 N36.25°, E100.52° 2984 JQ364719-JQ364721 
9 Heimahe, Qinghai (HH) YPX1018—1022 5 N36.51°, E99.55° 3181 JQ364722-JQ364726 
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Figure 1 Geographical sampling of P. raddei. Shading in bars corresponds to the major maternal lineages as shown in Figure 3. Locality 


abbreviations are indicated in Table 1. 


CytbA-L (5'-GAATY GGRGG WCAAC CAGTA 
GAAGA CCC-3') and Control B-H (5'-G TCCA TTGGA 
GGTTA AGATC TACCA-3') well amplified a 1 183 bp 
fragment of the partial sequence. Reactions were carried 
out in volumes of 25 ul as follows: 3 ul of 10 x buffer, 
17.875 ul dd H,O, 1 ul dNTP, 1 ul BSA, 0.125 ul Taq 
DNA polymerase, 1 ul template, and 1 ul each primer. 
PCR amplifications were run on PTC-200 thermal cycler 
(MJ Research Inc., USA) according to the following 
protocol: 94°C pre-denaturation for 5 min; 32 cycles of 
denaturation at 94°C for 1 min, annealing at 60°C for 1 
min, extension at 72°C for 1 min, and a final extension 
at 72°C for 10 min. After purification by a Clean-up Kit 
(V-gene), PCR products were sequenced on an ABI 3700 
automated sequencer, producing pherograms that were 
consistently readable. A 1 183 bp fragment was obtained 
for all specimens. 

A total of 185 new control region sequences of P. 
raddei plus the outgroup taxa were aligned using ClustalX 
(Thompson et al., 1997). We calculated the statistics Fu’s 
(1997) D* and F*, and Tajima’s D (Tajima, 1989) using 
DnaSP 4.10.8 (Rozas et al., 1999, 2003) to test the null 
hypothesis that all mutations in P. raddei were selectively 
neutral (Kimura, 1983). 


2.3 Matrilineal genealogy Phylogenetic and population 
genetic analyses were applied to investigate the 
matrilineal history and variation of P. raddei. Phylogenetic 
analysis of control region sequences including five 
outgroups was performed to reconstruct bifurcating 
trees using neighbor-joining (NJ) approaches using 


PAUP* v4b10 (Swofford, 2002). An unrooted Bayesian 
inference (BI) trees were constructed using MrBayes v3.1 
(Huelsenbeck and Ronquist, 2001). Modeltest 3.7 (Posada 
and Crandall, 1998) was used to select the optimal model 
of nucleotide evolution assuming the BIC criterion. BI 
analyses used one cold and three heated chains. The 
model employed six substitution types (Nst = 6) with 
a proportion of sites being invariant (rates = gamma). 
Samples of trees and parameters were drawn every | 000 
steps from a total of 10 000 000 MCMC generations. The 
first 2 000 000 generations (2 000 trees) were discarded as 
burn-in (Chain had not become stationary) and only the 
results from the last 8 000 000 generations (8 000 trees) 
were used to compute a majority rule consensus tree. 
The frequency of nodal resolution was termed a posterior 
probability. We ran three additional analyses starting with 
random trees. The final 32 000 post-burn-in trees from all 
four runs were used to construct the final majority rule 
tree. 


2.4 Estimates of genetic structure and historical 
demography Geographic structure for genealogical 
relationships was examined by mapping the lineages 
with locations. This determined the extent of geographic 
concordance among lineages. 

The significance of the correlation between genetic 
distance and geographic distance (isolation-by-distance) 
was used for the Mantel test (Mantel, 1967), and 
implemented in IBD 1.52 (Bohonak, 2002) with 10 000 
matrix randomizations. Hierarchical analyses of molecular 
variance (AMOVA; Excoffier et al., 1992) conducted 
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using Arlequin 3.0 (Excoffier et al., 2005) were employed 
to characterize population structure and genetic variation. 
Permutation tests of significance with 10 000 random 
permutations were used to determine genetic variance 
against null distributions. For the hierarchical analysis, 
populations were grouped according to geographic 
proximity. Nucleotide diversity, haplotype diversity, and 
pairwise measures of F,-values for all populations and 
each lineage were also calculated. 

Demographic patterns were examined using mismatch 
distributions. The observed vs. expected distributions 
of pairwise nucleotide differences between D-loop 
haplotypes were used to evaluate hypotheses of recent 
population growth. The fit of the observed data to the 
expected data (under the sudden growth model) was 
compared using the sum of squares deviations (SSD) 
between the observed and expected data estimated from 
2 000 parametric bootstrap replicates. Harpending’s 
raggedness index (HRI; Harpending, 1994) and its 
significance were calculated to quantify the smoothness of 
the observed mismatch distribution. Tajima’s D and Fu’s 
Fs tests were used to further test for possible demographic 
expansions. Unimodal mismatch distributions and 
negative values of Fs and D were assumed to indicate 
recent range expansions following recovery from a 
bottleneck (Tajima, 1989; Slatkin and Hudson, 1991; 
Rogers and Harpending, 1992; Fu, 1997). 

To estimate the shape of population growth of P. 
raddei through time, we constructed Bayesian skyline 
plots (BSP) implemented in BEAST v1.52 (Drummond et 
al., 2005; Drummond and Rambaut, 2006). We performed 
five Markov chain Monte Carlo (MCMC) runs for 50 
million iterations while sampling the genealogy and 
population size parameters every 2 000 iterations and 
discarding the first 10% as burn-in. The GTR + G+ I 
nucleotide substitution model, the most flexible model 
available for each lineage in BEAST, was used to allow 
the sample space of the parameters to be fully exploited 
(Huelsenbeck and Rannala, 2004). Although the mean 
substitution rate was fixed by assuming a conventional 
molecular clock for toads (Stéck et al., 2006), we used 
uncorrelated lognormal model (Drummond et al., 2006) 
to account for rate variation among lineages. Default 
settings of Bayesian priors were used. Demographic 
history through time was reconstructed using Tracer 1.5 
(Rambaut and Drummond, 2003). Time and effective 
population size were defined as millions of years and 0 
(NeT; T, generation time) for the Bayesian skyline plot, 
although this method cannot precisely estimate effective 
population size. 


2.5 Hypothesis testing In order to infer current and 
historical demographics of P. raddei, we implemented 
Hey and Nielsen’s (2004) isolation with migration model 
(Nielsen and Wakeley, 2001) using the application 
IMa (Hey and Nielsen, 2007). This method allowed us 
to simultaneously estimate the following parameters: 
effective female population size (N., = 9,/2 ug; where, x = 
W or E for the western and eastern regions, respectively, 
and ug = per-generation mutation rate) for each region (Ow 
and 6, are populations west and east of the mitochondrial 
discontinuity, respectively, as illustrated in Figure 1); 
ancestral population size (04) before population splitting; 
number of migrant females per generation (e. g., Myo. = 
Owmy/2); and time, in years, since population divergence 
(t = t/u; where t = time parameter in terms of u, and p is 
mutation rate per locus per year). For these calculations 
we assumed a per-generation mutation rate (uo) of 3.83 x 
10 ° (Stéck et al., 2006; Zhou and Song, 1997; generation 
time = 3 years). 

Five initial runs of IMa were used to obtain broad 
prior estimates of maximum parameter values that were 
specified in subsequent runs. All runs implemented 
Metropolis coupling over three chains with a single 
genealogy update per step. Analyses included a burn-in 
period of 500 000 steps. We ran four separate analyses 
with identical conditions, while using different random 
number as seeds. Runs were considered complete when 
effective sample sizes (ESS) for each parameter exceeded 
500. 

Priori, geologic and climatic processes predict two 
population histories: 1) a star-like population tree was 
assumed to reflect a fragmented large continuous range 
(Knowles, 2001). This hypothesis predicts that the 
divergence time of populations should date at most to 
the last glacial maxima (Figure 2 A); and 2) regional 
groupings reflected either the invasion from multiple 
refugia or long-term geographical isolation among 
regions (Figure 2 B). Fit of the D-loop gene tree to these 
two hypotheses was tested using the program Mesquite 
(Knowles and Maddison, 2002; Maddison and Maddison, 
2008). The S-statistic of Slatkin and Maddison (1989) 
provided a measure of fit between the gene tree and these 
alternative population trees. 1 000 coalescent genealogies 
were generated under each historical scenario and the 
distribution of S, the minimum number of sorting events 
required to explain the population subdivision (Slatkin 
and Maddison 1989), was recorded. The S value of our 
Bayes inferred genealogy was compared with the S values 
of the simulated genealogies to evaluate model fit. The 
test S-statistic was considered significant if it was found 
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to be less than 5% of values generated at random. 
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Figure 2 Population trees used in coalescent simulation tests of 
incomplete lineage sorting using Mesquite. Branch widths are 
indicated by specific @ values that are discussed in the text. Terminal 
branch widths sum to their representative Os (e. g., 0\/Py indicates 
that the width of each branch is the total and for the W group divided 
by the number of W populations; 0, is the ancestral population size), 
and branch lengths sum to ¢ (in generations). 


3. Results 


3.1. Sequence variation Control region sequences 
of 1 183 nucleotide positions in length were obtained 
from 185 ingroup individuals of P. raddei. Nucleotide 
compositions in P. raddei were as follows: C = 22.7%; T 
= 30.7%; A= 36.9%; and G = 9.7%. Of the positions, 162 
sites (13.7%) were variable, of which 108 (9.1%) were 
potentially phylogenetically informative. The transition/ 
transversion ratio was 2.6. Haplotype diversity (h) was 
0.99 and nucleotide diversity (pi) was 0.015. The negative 
value of Tajima’s D was not significant (-1.118, P > 0.10), 
showing no evidence that selection acted on this region. A 
neutral model of evolution could not be rejected. Within 
P. raddei, a total of 116 haplotypes were defined from all 
ingroup individuals. Four haplotypes (3.4%) were shared 
by some individuals from different sampling localities, 25 
(21.6%) were shared within the same sampling localities, 
and 87 (75.0%) were unique. A pattern of isolation-by- 
distance (IBD) was revealed by the Mantel test (r = 0.309, 
P = 0.026). The pairwise P distance ranged from 0.000 to 
0.035, with an average value being 0.015. 


3.2 Matrilineal relationships Given a low ratio of 
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potentially informative characters to taxa (108/185), 
many unresolved nodes were expected to be resolved. 
The BI tree resolved two major lineages, A and B, which 
were also resolved in the NJ tree (Figure 3). Lineage A 
(West lineage, 100% posterior probability) was confined 
to a small area in northwestern China. Lineage B, which 
had a poster probability of 100%, included the remaining 
haplotypes. This lineage was subdivided further into at 
least two sublineages. The toads from western and central 
(Midwest) China formed Sublineage B1 (Central lineage), 
and this lineage had a posterior probability of 100%. 
Sublineage B2 (East lineage) comprised individuals from 
northeastern China. The distribution of populations from 
the two main lineages overlapped in northwestern China. 
The hyplotypes from localities JY, YZ, and GY (All 
the abbreviations for localities here and after are fully 
provided in Table 1) clusted in both West and Central 
lineages, whereas the individuals from BX existed in all 
three lineages (Figures 1, 3). Differed from the NJ tree in 
only one way that some individuals from East populations 
(BL, HB, HL, XL, JN) branched by itself from the 
unresolved central node rather than grouping at the clade 
clustered by other individuals from East populations 
(Figure 3 B). 


3.3 Historical demography Mismatch frequencies were 
calculated separately for the three major lineages and 
each had a unimodal distribution. SSD and HRI indicated 
that the distributions for the eastern lineage were not 
significantly different from the distributions expected 
under the model of population expansion. Fu’s F and 
Tajima’s D statistical tests were negative for all three 
lineages, but the values for eastern lineage only were 
significant (P < 0.05). 

Estimated BSP for each lineage indicated that 10 
million generations were adequate. The BSP (Figure 4) 
resolved similar logistic population growth curves among 
lineages prior to 0.1 Ma and detected slightly recent 
population declines for the central lineages. During the 
Mid-late Pleistocene (0.4—0.5 Ma), the eastern lineage 
began to increase rapidly and then plateaued at about 0.2 
Ma. The size remained relatively stable during the early 
last glaciation until around 0.1 Ma when another rapid 
growth phase ensued. 


3.4 Coalescent analysis Demographic parameters 
indicated that the matrilineal effective population size 
of E [90% highest posterior density interval (HPD) 
1.54-1.77x10°] was similar to that of W (90% HPD 
1.41-1.78x10°). The time of isolation between these 
two regions was estimated to be 0.26 million years ago 
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Figure 3 Trees showing overall similarity and hypothesis of matrilineal relationships of P. raddei. A: Phylogenetic relationships of all 
haplotypes based on the D-loop sequences derived from NJ analysis. Numbers are NJ bootstraps values, only showing those higher than 
50%; B: Bayesian 50% majority-rule consensus network for 185 samples of P. raddei assuming the root resolved in the NJ tree. Posterior 
probabilities based on 10 000 post burn-in trees are shown above the branches; Letters or numbers are abbreviations for localities and 
individuals (See Table 1 for details); Major lineages are indicated by bars and corresponding names. Shading in bars corresponds to the 


geographic distributions in Figure 1. 
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(Ma) during the Mid-late Pleistocene and long before 
the last glacial maximum (LGM). Results from multiple, 
independent runs converged on values within this range, 
although HPD intervals were large and the distributions 
did not approach 0. Migration rates suggest that gene flow 
was mainly from west to east (My_,, = 4.2; Mg. w = 1.8). 

Either a large effective population size with high levels 
of gene flow or, more commonly, secondary contact 
between divergent lineages possibly has been attributed 
to a broad sympatry of evolutionarily divergent lineages. 
Consequently, we tested the null hypothesis of historical 
fragmentation and recent secondary contact. Rejection 
required discovery of a genetic signature of contiguous 
expansion consistent with one of two hypotheses of either 
incomplete lineage sorting resulting from fragmentation 
of a widespread ancestor or panmixia, alternatively. 
Consequently, we specified a population tree that reflected 
the fragmentation of a widespread ancestor (Figure 
2 A) in which it was specified as either 0 (panmixia) 
or as equal to the empirical estimate of divergence 
time obtained from the isolation with migration (IMa) 
analysis. Mesquite’s tree-search algorithm (Maddison 
and Maddison, 2004) was used to simulate gene trees by 
neutral coalescence within this population tree with an 
empirical estimate of divergence times (as measured in 
generations) and population sizes (Ne; both obtained from 
IMa). Simulated gene trees were then fit to the alternative 
population model of refugial isolation (Figure 2 B). 

All simulated gene trees within a population tree 
had an ancestral population size of N, = 6.8*10° (0, = 
26.04) applied to the root node, and we assigned E- and 
W-specific population sizes based on empirical estimates 
obtained from the IMa analysis, and those totals were 
divided equally among the respective terminals. For the 
recovered gene tree, Slatkin and Maddison’s s = 92. The 
discord predicted by coalescent simulation (average s 
= 89.3) could not reject the null hypothesis of refugial 
isolation. Conversely, the analysis rejected the alternative 
hypotheses for fragmentation either a widespread ancestor 
or panmixia (P = 0.02). 


4. Discussion 


4.1 Matrilineal structure and multiple refugia for P. 
raddei Excluding anthropogenic influences, either IBD 
and/or geography may have played an important role 
in the evolutionary history of amphibians and reptiles 
with limited vagility (Kuchta and Tan, 2005). IBD 
appears to be a factor involved in the population genetic 
structuring of P. raddei. Our genealogical analyses reveal 
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two lineages within the Chinese range of P. raddei: the 
eastern (A) and the western lineages (B) (Figure 3). 
Geographically, Lineage A covers most of the sampling 
region and Lineage B is restricted to the northeastern 
margin of the Qinghai-Tibet Plateau, the southwestern 
distribution of the species’ range (Figure 1). High level of 
mitochondrial divergence (uncorrected average divergence 
between the two major clades: 2.64 %; Nei and Tajima, 
1983), such as occurs in P. raddei, is often explained 
as the genetic signature of secondary contact between 
lineages previously isolated by the Pleistocene glaciation 
(Zhang et al., 2008; Ding et al., 2011). Based on the last 
connection of Africa and Europe and the fossil record, 
Stéck et al. (2006) suggest a divergence rate of 1.278%/Ma 
for D-loop within anurans. When this divergence rate 
is applied to our data, we date the initial split to about 
2.1 Ma, similar to that of another poikilotherm isolation 
event (Ding et al., 2011). This dating must be taken with 
caution because of the lack of calibration for P. raddei 
and the potential overestimation of timing recent events 
(Ho and Larson, 2006). Regardless, the date provides an 
approximate time frame. This divergence time is roughly 
congruent with the beginning of the Quaternary (2.5-1.8 
Ma), and environmental changes during this period 
may have driven the allopatry. At times of maximum 
glaciation, the Tibetan Plateau is known to have not 
been covered by ice sheets, yet periglacial environments 
contained valley-pediment glaciers (Shi et al., 1995). 
Lineage B from the northeastern margin of the Qinghai- 
Tibet Plateau (localities HH and HZ, Table 1) appears 
to have been isolated from Lineage A by mountains 
covered with ice sheets. It appears to have survived in the 
lower valley during ice ages, and subsequently dispersed 
during interglacial times. Unimodal distributions with 
low raggedness values and significant Fs values strongly 
supported the possibility of a sudden expansion in the 
eastern lineage. 

In contrast to the deep genealogical partitions 
commonly found in European and North American 
species (Hewitt, 2000), noticeable divisions do not occur 
in the central range of P. raddei. The observation of broad 
sympatry of divergent lineages suggests complex histories 
involving both allopatric isolation among refugia and 
prominent patterns of dispersal. 

Genealogical history is often used to draw conclusions 
concerning population history. When incomplete lineage 
sorting occurs, which is likely to occur in recently isolated 
populations, population history becomes difficult to 
infer (Charlesworth et al., 2003). Coalescent simulations 
provide a means of inferring population history, even 
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Figure 4 Mismatch distributions (left) and Bayesian skyline plots (right) depicting the demographic history for each lineage of P. raddei 
based on sequences from the mtDNA D-loop gene. For the mismatch distributions, the abscissa shows the number of pairwise differences 
between the compared haplotypes. The ordinate shows the proportion for each value. The bars show the observed frequency distribution, 
while the curve shows the distribution expected under the sudden-expansion model. For the skyline plots, data were analyzed for 20 million 
generations in the program BEAST v1.52 (Drummond and Rambaut, 2006; Drummond et al., 2005) using the appropriate model of evolution 
determined from ModelTest 3.7 (Posada and Crandall, 1998). The central line represents the median value for the log10 of the population 
size (Ne * t) and the dashed lines represent the upper and lower 95% credible intervals. 


in the face of incomplete lineage sorting, by assessing 
stochasticity inherent in the coalescent process. Our 
coalescent simulations are consistent with the scenario 
of isolation in at least western and eastern refugia during 
the Pleistocene glaciations. One of the refugia is located 
in the eastern monsoon region and the other one in lower 
elevations of the northeastern margin of the Qinghai- 
Tibet Plateau. Secondary contact between two previously 
isolated lineages appears to have occurred. This mixture 


of lineages is reflected in the AMOVA results, which 
indicate that the 47.3% of the observed genetic variation 
is due to and within population variation. 


4.2 Historical demography Assuming gene flow is 
freely occurring between the two major lineages, the 
focus turns to the demographic history of each maternal 
lineage. The unimodal mismatch distributions with low 
raggedness values and significantly negative Fu’s Fs 


No. 2 


strongly support the possibility of a recent expansion 
by the eastern lineage. Although D, Fs, and mismatch 
distributions provide insights into whether or not 
population growth has been expansive, they cannot 
discern the shape of population growth over time. Non- 
significant negative values of D and Fs neither reject the 
null hypothesis of population stability nor support the 
alternative hypothesis of expansive growth. However, 
such values cannot identify whether populations are 
expanding slowly, contracting, or remaining relatively 
constant in size. 

Bayesian skyline plots can provide insights into 
matrilineal patterns. This technique estimates N, through 
time and it does not require a specified demographic 
model (e. g., constant size, exponential growth, logistic 
growth, or expansive growth) prior to the analysis. 
Although the molecular clock for the toad may not be 
accurate, the historical population trend inferred from the 
Bayesian skyline plots seems similar to those of the other 
Chinese taxa since the Mid-late Pleistocene (Figure 4). 
During the last 100 000-150 000 years, populations of P. 
raddei from all lineages seem to have increased rapidly, 
and then plateaued at about 20 000 years ago, especially 
in the eastern and western lineages. Moreover, it seems 
that expansions occurred twice in the eastern group and 
probably much earlier than in the western one, suggesting 
more complicated changes of climates in the eastern part 
of northern China. 

Unlike species in temperate North America and 
Europe, which commonly show expansions in the post- 
LGM era (Mila et al., 2000, 2006, 2007; Kvist et al., 
2004; Hull and Girman, 2005; Fontanella et al., 2008; 
Guiher and Burbrink, 2008), population growth for 
lineages of P. raddei seems to have occurred prior to 
0.1-0.2 Ma. There is no evidence of population declines 
since the Late Pleistocene (Figure 4). Thus, populations 
do not appear to have been adversely affected by the Late 
Pleistocene to Holocene glaciations. The development of 
monsoons in East Asia since the Mid-late Pleistocene (Liu 
and Li, 1996; Zhang, 1999, 2002) is likely responsible for 
increasing available habitat for P. raddei. Paleoclimatic 
reconstructions based on d18 O values and pollen data 
show that the most extensive glaciation occurred during 
the Marine Isotope Stages 16-18 (MISI6—MIS 18; 0.6- 
0.7 Ma) in China (Wu et al., 2002). Subsequent climate 
oscillations seem to have been moderate in eastern China, 
where the populations of P. raddei were stably growing 
throughout the LGM. Recent studies show similar 
demographic patterns in eastern Chinese non-avian 
reptiles, frogs and birds; population expansions seem to 
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have begun prior to the LGM (Huang et al., 2007; Zhang 
et al., 2008; Li et al., 2009; Song et al., 2009; Ding et al., 
2011). These correlations further support our conclusions. 
In addition, the adaptation (or exaptation) of P. raddei to 
cold temperatures and desert habitats (Zhao and Adler, 
1993) might have permitted the rapid colonization of new 
areas and promoted population expansions. 
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